// reconstructs material interfaces from a mesh containing mixed cells (where several materials are mixed)
// this implementation is based on the youngs algorithm, generalized to arbitrary cell types and works
// on both 2D and 3D meshes. the main advantage of the youngs algorithm is it guarantees the material volume correctness.
// for 2D meshes, the AxisSymetric flag allows to switch between a pure 2D (plannar) algorithm and an axis symetric 2D algorithm
// handling volumes of revolution.

// Contributed by Thierry Carrard from
// Commissariat a l'Energie Atomique, (CEA)
// BP12, 91297 Arpajon, France


#include "vtkYoungsMaterialInterface.h"

#include <vtkCell.h>
#include <vtkEmptyCell.h>
#include <vtkPolygon.h>
#include <vtkConvexPointSet.h>
#include <vtkDataSet.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <vtkDataArray.h>
#include <vtkUnsignedCharArray.h>
#include <vtkIdTypeArray.h>
#include <vtkPoints.h>
#include <vtkCellArray.h>
#include <vtkUnstructuredGrid.h>
#include <vtkRectilinearGrid.h>
#include <vtkDoubleArray.h>
#include <vtkFloatArray.h>
#include <vtkMath.h>
#include <vtkIdList.h>

#include <vector>
#include <map>
#include <algorithm>

#include <assert.h>

#include <DebugStream.h>

#include "CellInterface.h"

#ifdef DEBUG
#define DBG_ASSERT(cond) assert(cond)
#else
#define DBG_ASSERT(cond) (void)0
#endif

// standard constructors and factory

//  Modifications:
//   Jeremy Meredith, Tue Dec  1 12:34:33 EST 2009
//   Made viable volume fraction range 2 orders of magnitude more exact.
//
vtkYoungsMaterialInterface::vtkYoungsMaterialInterface()
{
    this->FillMaterial = false;
    this->InverseNormal = false;
    this->AxisSymetric = false;
    this->OnionPeel = false;
    this->ReverseMaterialOrder = false;
    this->UseFractionAsDistance = false;
    this->VolumeFractionRange[0] = 0.0001;
    this->VolumeFractionRange[1] = 0.9999;
    this->TwoMaterialsOptimization = false;
}

vtkYoungsMaterialInterface::~vtkYoungsMaterialInterface()
{
}

void vtkYoungsMaterialInterface::SetNumberOfMaterials(int n)
{
    this->Materials.resize(n);
}

int vtkYoungsMaterialInterface::GetNumberOfMaterials()
{
    return this->Materials.size();
}

void vtkYoungsMaterialInterface::AddMaterial( const char* volume, const char* normal, const char* ordering )
{
    MaterialDescription md;
    md.volume = volume;
    md.normal = normal;
    md.normalX = "";
    md.normalY = "";
    md.normalZ = "";
    md.ordering = ordering;
    this->Materials.push_back( md );
}

/*
void vtkYoungsMaterialInterface::AddMaterial( const char* volume, const char* normalX, const char* normalY, const char* normalZ, const char* ordering )
{
  MaterialDescription md;
  md.volume = volume;
  md.normal = "";
  md.normalX = normalX;
  md.normalY = normalY;
  md.normalZ = normalZ;
  md.ordering = ordering;
  this->Materials.push_back( md );
}
*/

void vtkYoungsMaterialInterface::RemoveAllMaterials()
{
    //this->Materials.clear();
    this->SetNumberOfMaterials(0);
}


// internal classes
struct vtkYoungsMaterialInterface_IndexedValue
{
    double value;
    int index;
    inline bool operator < ( const vtkYoungsMaterialInterface_IndexedValue& iv ) const { return value < iv.value; }
};

struct vtkYoungsMaterialInterface_Mat
{
    // input
    vtkDataArray* fractionArray;
    vtkDataArray* normalArray;
    vtkDataArray* normalXArray;
    vtkDataArray* normalYArray;
    vtkDataArray* normalZArray;
    vtkDataArray* orderingArray;

    // temporary
    vtkIdType numberOfCells;
    vtkIdType numberOfPoints;
    vtkIdType cellCount;
    vtkIdType cellArrayCount;
    vtkIdType pointCount;
    vtkIdType* pointMap;

    // output
    std::vector<unsigned char> cellTypes;
    std::vector<vtkIdType> cells;
    vtkDataArray** outCellArrays;
    vtkDataArray** outPointArrays; // last point array is point coords
};


static inline void
vtkYoungsMaterialInterface_GetPointData(
                                        int nPointData,
                                        vtkDataArray** inPointArrays,
                                        vtkDataSet * input,
                                        std::vector< std::pair<int,vtkIdType> > & prevPointsMap,
                                        int nmat,
                                        vtkYoungsMaterialInterface_Mat * Mats ,
                                        int a, vtkIdType i, double* t)
{
    if( (i) >= 0 )
    {
        if(a<(nPointData-1))
        {
            DBG_ASSERT( /*i>=0 &&*/ i<inPointArrays[a]->GetNumberOfTuples());
            inPointArrays[a]->GetTuple( i , t );
        }
        else
        {
            DBG_ASSERT( a == (nPointData-1) );
            DBG_ASSERT( /*i>=0 &&*/ i<input->GetNumberOfPoints());
            input->GetPoint( i , t );
        }
    }
    else
    {
        int j=-i-1;
        DBG_ASSERT(j>=0 && j<prevPointsMap.size());
        int prev_m = prevPointsMap[j].first;
        DBG_ASSERT(prev_m>=0 && prev_m<nmat);
        vtkIdType prev_i = ( prevPointsMap[j].second );
        DBG_ASSERT(prev_i>=0 && prev_i<Mats[prev_m].outPointArrays[a]->GetNumberOfTuples());
        Mats[prev_m].outPointArrays[a]->GetTuple( prev_i , t );
    }
}

#define GET_POINT_DATA(a,i,t) vtkYoungsMaterialInterface_GetPointData(nPointData,inPointArrays,input,prevPointsMap,nmat,Mats,a,i,t)

//  Modifications:
//   Jeremy Meredith, Thu Mar 25 11:30:20 EDT 2010
//   Renamed.  There was a namespace conflict and the compiler got confused
//   at runtime (but not compile time, frustratingly).
struct YoungsCellInfo
{
    double points[vtkYoungsMaterialInterface::MAX_CELL_POINTS][3];
    vtkIdType pointIds[vtkYoungsMaterialInterface::MAX_CELL_POINTS];
    int triangulation[vtkYoungsMaterialInterface::MAX_CELL_POINTS*4];
    int edges[vtkYoungsMaterialInterface::MAX_CELL_POINTS][2];

    int dim;
    int np;
    int nf;
    int ntri;
    int type;
    int nEdges;

    bool triangulationOk;
    bool needTriangulation;

    inline YoungsCellInfo() : dim(2), np(0), nf(0), ntri(0), type(VTK_EMPTY_CELL), nEdges(0), triangulationOk(false), needTriangulation(false) {}
};

int vtkYoungsMaterialInterface::CellProduceInterface( int dim, int np, double fraction, double minFrac, double maxFrac )
{
    return
        (
         (dim==3 && np>=4) ||
         (dim==2 && np>=3)
         ) &&
        (
         this->UseFractionAsDistance ||
         (
          ( fraction > minFrac ) &&
          ( fraction < maxFrac || this->FillMaterial )
          )
         ) ;

}

//
//  Modifications:
//    Kathleen Biagas, Thu Aug 11, 2022
//    Support VTK 9's change in storage for vtkCellArray class:
//    separate connectivity and offsets arrays.
//
int vtkYoungsMaterialInterface::Execute(vtkDataSet *input,
                                        vtkDataSet **output)
{
    // variables visible by debugger
    static const double minFrac = this->VolumeFractionRange[0];
    static const double maxFrac = this->VolumeFractionRange[1];

    int nmat = this->Materials.size();
    int nCellData = input->GetCellData()->GetNumberOfArrays();
    int nPointData = input->GetPointData()->GetNumberOfArrays();
    vtkIdType nCells = input->GetNumberOfCells();
    vtkIdType nPoints = input->GetNumberOfPoints();

    // -------------- temporary data initialization -------------------
    vtkDataArray** inCellArrays = new vtkDataArray* [nCellData] ;
    for(int i=0;i<nCellData;i++)
    {
        inCellArrays[i] = input->GetCellData()->GetArray(i);
    }

    vtkDataArray** inPointArrays = new vtkDataArray* [nPointData+1]; // last point array is point coords
    int* pointArrayOffset = new int [ nPointData+1 ];
    int pointDataComponents = 0;
    for(int i=0;i<nPointData;i++)
    {
        inPointArrays[i] = input->GetPointData()->GetArray(i);
        pointArrayOffset[i] = pointDataComponents;
        pointDataComponents += inPointArrays[i]->GetNumberOfComponents();
    }
    // we add another data array for point coords
    pointArrayOffset[nPointData] = pointDataComponents;
    pointDataComponents += 3;
    nPointData++;

    int origCoordType = VTK_FLOAT;
    if (input->GetDataObjectType() == VTK_RECTILINEAR_GRID)
    {
        origCoordType = vtkRectilinearGrid::SafeDownCast(input)->GetXCoordinates()->GetDataType();
    }
    else
    {
        origCoordType = vtkPointSet::SafeDownCast(input)->GetPoints()->GetDataType();
    }

    vtkYoungsMaterialInterface_Mat* Mats = new vtkYoungsMaterialInterface_Mat[nmat];
    {
        int m=0;
        for( std::vector<MaterialDescription>::iterator it=this->Materials.begin(); it!=this->Materials.end(); ++it , ++m )
        {
            Mats[m].fractionArray = input->GetCellData()->GetArray( (*it).volume.c_str() );
            Mats[m].normalArray = input->GetCellData()->GetArray( (*it).normal.c_str() );
            Mats[m].normalXArray = input->GetCellData()->GetArray( (*it).normalX.c_str() );
            Mats[m].normalYArray = input->GetCellData()->GetArray( (*it).normalY.c_str() );
            Mats[m].normalZArray = input->GetCellData()->GetArray( (*it).normalZ.c_str() );
            Mats[m].orderingArray = input->GetCellData()->GetArray( (*it).ordering.c_str() );

            if( Mats[m].fractionArray == 0 )
            {
                debug3<<"Material "<<m<<": volume fraction array '"<<(*it).volume<<"' not found\n";
            }
            if( Mats[m].orderingArray == 0 )
            {
                debug3<<"Material "<<m<<" material ordering array '"<<(*it).ordering<<"' not found\n";
            }
            if( Mats[m].normalArray==0 && Mats[m].normalXArray==0 && Mats[m].normalYArray==0 && Mats[m].normalZArray==0 )
            {
                debug3<<"Material "<<m<<" normal  array '"<<(*it).normal<<"' not found\n";
            }


            Mats[m].numberOfCells = 0;
            Mats[m].cellCount = 0;
            Mats[m].cellArrayCount = 0;

            Mats[m].outCellArrays = new vtkDataArray* [ nCellData ];
            for(int i=0;i<nCellData;i++)
            {
                Mats[m].outCellArrays[i] = vtkDataArray::CreateDataArray( inCellArrays[i]->GetDataType() );
                Mats[m].outCellArrays[i]->SetName( inCellArrays[i]->GetName() );
                Mats[m].outCellArrays[i]->SetNumberOfComponents( inCellArrays[i]->GetNumberOfComponents() );
            }

            Mats[m].numberOfPoints = 0;
            Mats[m].pointCount = 0;
            Mats[m].outPointArrays = new vtkDataArray* [ nPointData ];

            for(int i=0;i<(nPointData-1);i++)
            {
                Mats[m].outPointArrays[i] = vtkDataArray::CreateDataArray( inPointArrays[i]->GetDataType() );
                Mats[m].outPointArrays[i]->SetName( inPointArrays[i]->GetName() );
                Mats[m].outPointArrays[i]->SetNumberOfComponents( inPointArrays[i]->GetNumberOfComponents() );
            }
            // JSM August  3, 2009: this was originally a double array, not float
            //Mats[m].outPointArrays[nPointData-1] = vtkFloatArray::New();
            Mats[m].outPointArrays[nPointData-1] = vtkDoubleArray::New();
            Mats[m].outPointArrays[nPointData-1]->SetName("Points");
            Mats[m].outPointArrays[nPointData-1]->SetNumberOfComponents(3);
        }
    }


    // --------------- per material number of interfaces estimation ------------
    for(vtkIdType c=0;c<nCells;c++)
    {
        vtkCell* vtkcell = input->GetCell(c);
        int cellDim = vtkcell->GetCellDimension();
        int np = vtkcell->GetNumberOfPoints();
        int nf = vtkcell->GetNumberOfFaces();

        for(int m=0;m<nmat;m++)
        {
            double fraction = ( Mats[m].fractionArray != 0 ) ? Mats[m].fractionArray->GetTuple1(c) : 0 ;
            if( this->CellProduceInterface(cellDim,np,fraction,minFrac,maxFrac) )
            {
                if( cellDim == 2 )
                {
                    Mats[m].numberOfPoints += 2;
                }
                else
                {
                    Mats[m].numberOfPoints += nf;
                }
                if( this->FillMaterial )
                {
                    Mats[m].numberOfPoints += np-1;
                }
                Mats[m].numberOfCells ++;
            }
        }
    }


    // allocation of output arrays
    for(int m=0;m<nmat;m++)
    {
        debug3<<"Mat #"<<m<<" : cells="<<Mats[m].numberOfCells<<", points="<<Mats[m].numberOfPoints<<", FillMaterial="<<this->FillMaterial<<"\n";
        for(int i=0;i<nCellData;i++)
        {
            Mats[m].outCellArrays[i]->Allocate( Mats[m].numberOfCells * Mats[m].outCellArrays[i]->GetNumberOfComponents() );
        }
        for(int i=0;i<nPointData;i++)
        {
            Mats[m].outPointArrays[i]->Allocate( Mats[m].numberOfPoints * Mats[m].outPointArrays[i]->GetNumberOfComponents() );
        }
        Mats[m].cellTypes.reserve( Mats[m].numberOfCells );
        Mats[m].cells.reserve( Mats[m].numberOfCells + Mats[m].numberOfPoints );
        Mats[m].pointMap = new vtkIdType[ nPoints ];
        for(vtkIdType i=0;i<nPoints;i++) { Mats[m].pointMap[i] = -1; }
    }

    // --------------------------- core computation --------------------------
    vtkIdList *ptIds = vtkIdList::New();
    vtkPoints *pts = vtkPoints::New(origCoordType);
    vtkConvexPointSet* cpsCell = vtkConvexPointSet::New();

    double* interpolatedValues = new double[ MAX_CELL_POINTS * pointDataComponents ];
    vtkYoungsMaterialInterface_IndexedValue * matOrdering = new vtkYoungsMaterialInterface_IndexedValue[nmat];

    std::vector< std::pair<int,vtkIdType> > prevPointsMap;
    prevPointsMap.reserve( MAX_CELL_POINTS*nmat );

    for(vtkIdType ci=0;ci<nCells;ci++)
    {
        int interfaceEdges[MAX_CELL_POINTS*2];
        double interfaceWeights[MAX_CELL_POINTS];
        int nInterfaceEdges;

        int insidePointIds[MAX_CELL_POINTS];
        int nInsidePoints;

        int outsidePointIds[MAX_CELL_POINTS];
        int nOutsidePoints;

        int outCellPointIds[MAX_CELL_POINTS];
        int nOutCellPoints;

        double referenceVolume = 1.0;
        double normal[3];

        prevPointsMap.clear();

        // sort materials
        int nEffectiveMat = 0;
        for(int mi=0;mi<nmat;mi++)
        {
            matOrdering[mi].index = mi;
            matOrdering[mi].value = ( Mats[mi].orderingArray != 0 ) ? Mats[mi].orderingArray->GetTuple1(ci) : 0.0;

            double fraction = ( Mats[mi].fractionArray != 0 ) ? Mats[mi].fractionArray->GetTuple1(ci) : 0 ;
            if( fraction>minFrac ) nEffectiveMat++;
        }
        std::stable_sort( matOrdering , matOrdering+nmat );

        bool twoMaterialOptimization =  this->TwoMaterialsOptimization && (nEffectiveMat==2);
        if( twoMaterialOptimization )
        {
            debug3<<"2 material optimization triggered for cell #"<<ci<<"\n";
        }

        // read cell information for the first iteration
        // a temporary cell will then be generated after each iteration for the next one.
        vtkCell* vtkcell = input->GetCell(ci);
        YoungsCellInfo cell;
        cell.dim = vtkcell->GetCellDimension();
        cell.np = vtkcell->GetNumberOfPoints();
        cell.nf = vtkcell->GetNumberOfFaces();
        cell.type = vtkcell->GetCellType();

        // copy points and point ids to lacal arrays.
        // IMPORTANT NOTE : A negative point id refers to a point in the previous material.
        // the material number and real point id can be found through the prevPointsMap.
        if (cell.type == VTK_PIXEL)
        {
            // JSM December 23, 2009:
            // Add PIXEL support; assumptions are made that the cell
            // connectivity is ordered clockwise like a POLYGON or QUAD.
            int s0 = 0, s1 = 1, s2 = 3, s3 = 2;
            cell.pointIds[0] = vtkcell->GetPointId(s0);
            vtkcell->GetPoints()->GetPoint( s0 , cell.points[0] );
            cell.pointIds[1] = vtkcell->GetPointId(s1);
            vtkcell->GetPoints()->GetPoint( s1 , cell.points[1] );
            cell.pointIds[2] = vtkcell->GetPointId(s2);
            vtkcell->GetPoints()->GetPoint( s2 , cell.points[2] );
            cell.pointIds[3] = vtkcell->GetPointId(s3);
            vtkcell->GetPoints()->GetPoint( s3 , cell.points[3] );
            cell.type = VTK_POLYGON;
        }
        else
        {
            for(int p=0;p<cell.np;p++)
            {
                cell.pointIds[p] = vtkcell->GetPointId(p);
                DBG_ASSERT( cell.pointIds[p]>=0 && cell.pointIds[p]<nPoints );
                vtkcell->GetPoints()->GetPoint( p , cell.points[p] );
            }
        }


        // Triangulate cell.
        // IMPORTANT NOTE: triangulation is given with mesh point ids (not local cell ids)
        // and are translated to cell local point ids.
        cell.needTriangulation = false;
        cell.triangulationOk = vtkcell->Triangulate(ci,ptIds,pts);
        cell.ntri = 0;
        if( cell.triangulationOk )
        {
            cell.ntri = ptIds->GetNumberOfIds() / (cell.dim+1);
            for(int i=0;i<(cell.ntri*(cell.dim+1));i++)
            {
                vtkIdType j = std::find( cell.pointIds , cell.pointIds+cell.np , ptIds->GetId(i) ) - cell.pointIds;
                DBG_ASSERT( j>=0 && j<cell.np );
                cell.triangulation[i] = j;
            }
        }
        else
        {
            debug3<<"Triangulation failed on primary cell\n";
        }

        // get 3D cell edges.
        if( cell.dim == 3 )
        {
            vtkCell3D* cell3D = vtkCell3D::SafeDownCast( vtkcell );
            cell.nEdges = vtkcell->GetNumberOfEdges();
            for(int i=0;i<cell.nEdges;i++)
            {
                const vtkIdType *edgePoints = nullptr;
                cell3D->GetEdgePoints(i,edgePoints);
                cell.edges[i][0] = edgePoints[0];
                DBG_ASSERT( cell.edges[i][0]>=0 && cell.edges[i][0]<cell.np );
                cell.edges[i][1] = edgePoints[1];
                DBG_ASSERT( cell.edges[i][1]>=0 && cell.edges[i][1]<cell.np );
            }
        }

        // for debugging : ensures that we don't read anything from cell, but only from previously filled arrays
        vtkcell = 0;

        int processedEffectiveMat = 0;

        // Loop for each material. Current cell is iteratively cut.
        for(int mi=0;mi<nmat;mi++)
        {
            int m = this->ReverseMaterialOrder ? matOrdering[nmat-1-mi].index : matOrdering[mi].index;

            // get volume fraction and interface plane normal from input arrays
            double fraction = ( Mats[m].fractionArray != 0 ) ? Mats[m].fractionArray->GetTuple1(ci) : 0 ;
            if (!twoMaterialOptimization)
            {
                fraction = (referenceVolume>0) ? (fraction/referenceVolume) : 0.0 ;
            }

            if( this->CellProduceInterface(cell.dim,cell.np,fraction,minFrac,maxFrac) )
            {
                YoungsCellInfo nextCell; // empty cell by default
                int interfaceCellType = VTK_EMPTY_CELL;

                if( mi==0 || ( !this->OnionPeel && !twoMaterialOptimization) )
                {
                    normal[0]=0; normal[1]=0; normal[2]=0;

                    if( Mats[m].normalArray != 0 ) Mats[m].normalArray->GetTuple(ci,normal);
                    if( Mats[m].normalXArray != 0 ) normal[0] = Mats[m].normalXArray->GetTuple1(ci);
                    if( Mats[m].normalYArray != 0 ) normal[1] = Mats[m].normalYArray->GetTuple1(ci);
                    if( Mats[m].normalZArray != 0 ) normal[2] = Mats[m].normalZArray->GetTuple1(ci);

                    // work-around for degenerated normals
                    if( vtkMath::Norm(normal) == 0.0 ) // should it be <EPSILON ?
                    {
                        //debug3<<"Nul normal\n";
                        normal[0]=1.0;
                        normal[1]=0.0;
                        normal[2]=0.0;
                    }
                    else
                    {
                        vtkMath::Normalize( normal );
                    }
                    if( this->InverseNormal )
                    {
                        normal[0] = -normal[0];
                        normal[1] = -normal[1];
                        normal[2] = -normal[2];
                    }
                }

                // 2 material case optimization
                if( twoMaterialOptimization && processedEffectiveMat>0 )
                {
                    normal[0] = -normal[0];
                    normal[1] = -normal[1];
                    normal[2] = -normal[2];
                    if( fraction > minFrac )
                    {
                        processedEffectiveMat ++;
                    }
                }

                // -= case where the entire input cell is passed through =-
                if( fraction>this->VolumeFractionRange[1] && this->FillMaterial )
                {
                    interfaceCellType = cell.type;
                    //Mats[m].cellTypes.push_back( cell.type );
                    nOutCellPoints = nInsidePoints = cell.np;
                    nInterfaceEdges = 0;
                    nOutsidePoints = 0;
                    for(int p=0;p<cell.np;p++) { outCellPointIds[p] = insidePointIds[p] = p; }
                    // remaining volume is an empty cell (nextCell is left as is)
                }

                // -= case where the entire cell is ignored =-
                else if( fraction<this->VolumeFractionRange[0] || (fraction>this->VolumeFractionRange[1] && !this->FillMaterial) || !cell.triangulationOk )
                {
                    interfaceCellType = VTK_EMPTY_CELL;
                    //Mats[m].cellTypes.push_back( VTK_EMPTY_CELL );

                    nOutCellPoints = 0;
                    nInterfaceEdges = 0;
                    nInsidePoints = 0;
                    nOutsidePoints = 0;

                    // remaining volume is the same cell
                    nextCell = cell;

                    if( !cell.triangulationOk )
                    {
                        debug3<<"Cell triangulation failed\n";
                    }
                }

                // -= 2D case =-
                else if( cell.dim == 2 )
                {
                    int nRemCellPoints;
                    int remCellPointIds[MAX_CELL_POINTS];

                    int triangles[MAX_CELL_POINTS][3];
                    for(int i=0;i<cell.ntri;i++) for(int j=0;j<3;j++)
                    {
                        triangles[i][j] = cell.triangulation[i*3+j];
                        DBG_ASSERT( triangles[i][j]>=0 && triangles[i][j]<cell.np );
                    }

                    bool interfaceFound = CellInterface::cellInterface2D(
                                                                         cell.points, cell.np,
                                                                         triangles, cell.ntri,
                                                                         fraction, normal,
                                                                         this->AxisSymetric,
                                                                         this->UseFractionAsDistance,
                                                                         interfaceEdges, interfaceWeights,
                                                                         nOutCellPoints, outCellPointIds,
                                                                         nRemCellPoints, remCellPointIds);

                    if( interfaceFound )
                    {
                        nInterfaceEdges = 2;
                        interfaceCellType = this->FillMaterial ? VTK_POLYGON : VTK_LINE;
                        //Mats[m].cellTypes.push_back( this->FillMaterial ? VTK_POLYGON : VTK_LINE );

                        // remaining volume is a polygon
                        nextCell.dim = 2;
                        nextCell.np = nRemCellPoints;
                        nextCell.nf = nRemCellPoints;
                        nextCell.type = VTK_POLYGON;

                        // build polygon triangulation for next iteration
                        nextCell.ntri = nextCell.np-2;
                        for(int i=0;i<nextCell.ntri;i++)
                        {
                            nextCell.triangulation[i*3+0] = 0;
                            nextCell.triangulation[i*3+1] = i+1;
                            nextCell.triangulation[i*3+2] = i+2;
                        }
                        nextCell.triangulationOk = true;
                        nextCell.needTriangulation = false;

                        // populate prevPointsMap and next iteration cell point ids
                        int ni = 0;
                        for(int i=0;i<nRemCellPoints;i++)
                        {
                            vtkIdType id = remCellPointIds[i];
                            if( id < 0 )
                            {
                                id = -(int)prevPointsMap.size()-1;
                                DBG_ASSERT( (-id-1) == prevPointsMap.size() );
                                prevPointsMap.push_back( std::make_pair( m , Mats[m].pointCount+ni ) ); // intersection points will be added first
                                ni++;
                            }
                            else
                            {
                                DBG_ASSERT( id>=0 && id<cell.np );
                                id = cell.pointIds[ id ];
                            }
                            nextCell.pointIds[i] = id;
                        }
                        DBG_ASSERT( ni == nInterfaceEdges );

                        // filter out points inside material volume
                        nInsidePoints = 0;
                        for(int i=0;i<nOutCellPoints;i++)
                        {
                            if( outCellPointIds[i] >= 0 ) insidePointIds[nInsidePoints++] = outCellPointIds[i];
                        }

                        if( ! this->FillMaterial ) // keep only interface points
                        {
                            int n = 0;
                            for(int i=0;i<nOutCellPoints;i++)
                            {
                                if( outCellPointIds[i] < 0 ) outCellPointIds[n++] = outCellPointIds[i];
                            }
                            nOutCellPoints = n;
                        }
                    }
                    else
                    {
                        debug3<<"no interface found for cell "<<ci<<", mi="<<mi<<", m="<<m<<", frac="<<fraction<<"\n";
                        nInterfaceEdges = 0;
                        nOutCellPoints = 0;
                        nInsidePoints = 0;
                        nOutsidePoints = 0;
                        interfaceCellType = VTK_EMPTY_CELL;
                        //Mats[m].cellTypes.push_back( VTK_EMPTY_CELL );
                        // remaining volume is the original cell left unmodified
                        nextCell = cell;
                    }

                }

                // -= 3D case =-
                else
                {
                    int tetras[MAX_CELL_POINTS][4];
                    for(int i=0;i<cell.ntri;i++) for(int j=0;j<4;j++)
                    {
                        tetras[i][j] = cell.triangulation[i*4+j];
                    }

                    // compute innterface polygon
                    CellInterface::cellInterface3D(
                                                   cell.np, cell.points,
                                                   cell.nEdges, cell.edges,
                                                   cell.ntri, tetras,
                                                   fraction, normal,
                                                   this->UseFractionAsDistance,
                                                   nInterfaceEdges, interfaceEdges, interfaceWeights,
                                                   nInsidePoints, insidePointIds,
                                                   nOutsidePoints, outsidePointIds );

                    if( nInterfaceEdges>cell.nf || nInterfaceEdges<3 ) // degenerated case, considered as null interface
                    {
                        nInterfaceEdges = 0;
                        nOutCellPoints = 0;
                        nInsidePoints = 0;
                        nOutsidePoints = 0;
                        interfaceCellType = VTK_EMPTY_CELL;
                        //Mats[m].cellTypes.push_back( VTK_EMPTY_CELL );

                        // in this case, next iteration cell is the same
                        nextCell = cell;
                    }
                    else
                    {
                        nOutCellPoints = 0;

                        for(int e=0;e<nInterfaceEdges;e++)
                        {
                            outCellPointIds[nOutCellPoints++] = -e -1;
                        }

                        if(this->FillMaterial)
                        {
                            interfaceCellType = VTK_CONVEX_POINT_SET;
                            //Mats[m].cellTypes.push_back( VTK_CONVEX_POINT_SET );
                            for(int p=0;p<nInsidePoints;p++)
                            {
                                outCellPointIds[nOutCellPoints++] = insidePointIds[p];
                            }
                        }
                        else
                        {
                            interfaceCellType = VTK_POLYGON;
                            //Mats[m].cellTypes.push_back( VTK_POLYGON );
                        }

                        // remaining volume is a convex point set
                        // IMPORTANT NOTE: next iteration cell cannot be entirely built right now.
                        // in this particular case we'll finish it at the end of the material loop
                        if( mi<(nmat-1) && !twoMaterialOptimization )
                        {
                            // JSM August  3, 2009: this next line was dropped between
                            // the two versions Thierry distributed.  It turns out to
                            // have been critical, causing bad 3D results if neglected.
                            nextCell.dim = 3;
                            nextCell.type = VTK_CONVEX_POINT_SET;
                            nextCell.np = nInterfaceEdges + nOutsidePoints;
                            vtkcell = cpsCell;
                            vtkcell->Points->Reset();
                            vtkcell->PointIds->Reset();
                            vtkcell->Points->SetNumberOfPoints( nextCell.np );
                            vtkcell->PointIds->SetNumberOfIds( nextCell.np );
                            for(int i=0;i<nextCell.np;i++)
                            {
                                vtkcell->PointIds->SetId( i, i );
                            }
                            // nf, ntri and triangulation have to be computed later on, when point coords are computed
                            nextCell.needTriangulation = true;
                        }

                        for(int i=0;i<nInterfaceEdges;i++)
                        {
                            vtkIdType id = -(int)prevPointsMap.size() - 1;
                            DBG_ASSERT( (-id-1) == prevPointsMap.size() );
                            prevPointsMap.push_back( std::make_pair( m , Mats[m].pointCount+i ) ); // we know that interpolated points will be added consecutively
                            nextCell.pointIds[i] = id;
                        }
                        for(int i=0;i<nOutsidePoints;i++)
                        {
                            nextCell.pointIds[nInterfaceEdges+i] = cell.pointIds[ outsidePointIds[i] ];
                        }
                    }

                    // check correctness of next cell's point ids
                    for(int i=0;i<nextCell.np;i++)
                    {
                        DBG_ASSERT( ( nextCell.pointIds[i]<0 && (-nextCell.pointIds[i]-1)<prevPointsMap.size() ) || ( nextCell.pointIds[i]>=0 && nextCell.pointIds[i]<nPoints ) );
                    }

                } // End 3D case


                //  create output cell
                if( interfaceCellType != VTK_EMPTY_CELL )
                {

                    // set type of cell
                    Mats[m].cellTypes.push_back( interfaceCellType );

                    // interpolate point values for cut edges
                    for(int e=0;e<nInterfaceEdges;e++)
                    {
                        double t = interfaceWeights[e];
                        for(int p=0;p<nPointData;p++)
                        {
                            double v0[16];
                            double v1[16];
                            int nc = Mats[m].outPointArrays[p]->GetNumberOfComponents();
                            int ep0 = cell.pointIds[ interfaceEdges[e*2+0] ];
                            int ep1 = cell.pointIds[ interfaceEdges[e*2+1] ];
                            GET_POINT_DATA( p , ep0 , v0 );
                            GET_POINT_DATA( p , ep1 , v1 );
                            for(int c=0;c<nc;c++)
                            {
                                interpolatedValues[ e*pointDataComponents + pointArrayOffset[p] + c ] = v0[c] + t * ( v1[c] - v0[c] );
                            }
                        }
                    }

                    // copy point values
                    for(int e=0;e<nInterfaceEdges;e++)
                    {
                        for(int a=0;a<nPointData;a++)
                        {
                            DBG_ASSERT( nptId == Mats[m].outPointArrays[a]->GetNumberOfTuples() );
                            Mats[m].outPointArrays[a]->InsertNextTuple( interpolatedValues + e*pointDataComponents + pointArrayOffset[a] );
                        }
                    }
                    int pointsCopied = 0;
                    int prevMatInterfToBeAdded = 0;
                    if( this->FillMaterial )
                    {
                        for(int p=0;p<nInsidePoints;p++)
                        {
                            vtkIdType ptId = cell.pointIds[ insidePointIds[p] ];
                            if( ptId>=0 )
                            {
                                if( Mats[m].pointMap[ptId] == -1 )
                                {
                                    vtkIdType nptId = Mats[m].pointCount + nInterfaceEdges + pointsCopied;
                                    Mats[m].pointMap[ptId] = nptId;
                                    pointsCopied++;
                                    for(int a=0;a<nPointData;a++)
                                    {
                                        DBG_ASSERT( nptId == Mats[m].outPointArrays[a]->GetNumberOfTuples() );
                                        double tuple[16];
                                        GET_POINT_DATA( a, ptId, tuple );
                                        Mats[m].outPointArrays[a]->InsertNextTuple( tuple );
                                    }
                                }
                            }
                            else
                            {
                                prevMatInterfToBeAdded++;
                            }
                        }
                    }

                    // populate connectivity array
                    // and add extra points from previous edge intersections that are used but not inserted yet
                    int prevMatInterfAdded = 0;
                    Mats[m].cells.push_back( nOutCellPoints ); Mats[m].cellArrayCount++;
                    for(int p=0;p<nOutCellPoints;p++)
                    {
                        int nptId;
                        int pointIndex = outCellPointIds[p];
                        if( pointIndex >=0 ) // an original point (not an edge intersection)
                        {
                            DBG_ASSERT( pointIndex>=0 && pointIndex<cell.np );
                            int ptId = cell.pointIds[ pointIndex ];
                            if( ptId>=0 )
                            {
                                DBG_ASSERT( ptId>=0 && ptId<nPoints ); // OUI, car interface d'une iteration precedente. que faire ...
                                nptId = Mats[m].pointMap[ptId];
                            }
                            else
                            {
                                nptId = Mats[m].pointCount + nInterfaceEdges + pointsCopied + prevMatInterfAdded;
                                prevMatInterfAdded++;
                                for(int a=0;a<nPointData;a++)
                                {
                                    DBG_ASSERT( nptId == Mats[m].outPointArrays[a]->GetNumberOfTuples() );
                                    double tuple[16];
                                    GET_POINT_DATA( a, ptId, tuple );
                                    Mats[m].outPointArrays[a]->InsertNextTuple( tuple );
                                }
                            }
                        }
                        else
                        {
                            int interfaceIndex = -pointIndex - 1;
                            DBG_ASSERT( interfaceIndex>=0 && interfaceIndex<nInterfaceEdges );
                            nptId = Mats[m].pointCount + interfaceIndex ;
                        }
                        DBG_ASSERT( nptId>=0 && nptId<(Mats[m].pointCount+nInterfaceEdges+pointsCopied+prevMatInterfToBeAdded) );
                        Mats[m].cells.push_back( nptId ); Mats[m].cellArrayCount++;
                    }

                    Mats[m].pointCount += nInterfaceEdges + pointsCopied + prevMatInterfAdded;

                    // copy cell arrays
                    for(int a=0;a<nCellData;a++)
                    {
                        Mats[m].outCellArrays[a]->InsertNextTuple( inCellArrays[a]->GetTuple(ci) );
                    }
                    Mats[m].cellCount ++;

                    // check for equivalence between counters and container sizes
                    DBG_ASSERT( Mats[m].cellCount == Mats[m].cellTypes.size() );
                    DBG_ASSERT( Mats[m].cellArrayCount == Mats[m].cells.size() );

                    // populate next iteration cell's point coords
                    for(int i=0;i<nextCell.np;i++)
                    {
                        DBG_ASSERT( ( nextCell.pointIds[i]<0 && (-nextCell.pointIds[i]-1)<prevPointsMap.size() ) || ( nextCell.pointIds[i]>=0 && nextCell.pointIds[i]<nPoints ) );
                        GET_POINT_DATA( (nPointData-1) , nextCell.pointIds[i] , nextCell.points[i] );
                    }

                    // for the convex point set, we need to first compute point coords before triangulation (no fixed topology)
                    if( nextCell.needTriangulation && mi<(nmat-1) )
                    {
                        // JSM August  3, 2009: this next loop was dropped between
                        // the two versions Thierry distributed.  It turns out to
                        // have been critical: leaving it out apparently causes UMRs
                        // and bad results, in particular slowdowns to 1 sec/cell.
                        for(int i=0; i<nextCell.np; i++)
                        {
                            vtkcell->Points->SetPoint(i, nextCell.points[i]);
                        }
                        vtkcell->Initialize();
                        nextCell.nf = vtkcell->GetNumberOfFaces();
                        if( nextCell.dim == 3 )
                        {
                            vtkCell3D* cell3D = vtkCell3D::SafeDownCast( vtkcell );
                            nextCell.nEdges = vtkcell->GetNumberOfEdges();
                            for(int i=0;i<nextCell.nEdges;i++)
                            {
                                const vtkIdType * edgePoints = nullptr;
                                cell3D->GetEdgePoints(i,edgePoints);
                                nextCell.edges[i][0] = edgePoints[0];
                                DBG_ASSERT( nextCell.edges[i][0]>=0 && nextCell.edges[i][0]<nextCell.np );
                                nextCell.edges[i][1] = edgePoints[1];
                                DBG_ASSERT( nextCell.edges[i][1]>=0 && nextCell.edges[i][1]<nextCell.np );
                            }
                        }
                        nextCell.triangulationOk = vtkcell->Triangulate(ci,ptIds,pts);
                        nextCell.ntri = 0;
                        if( nextCell.triangulationOk )
                        {
                            nextCell.ntri = ptIds->GetNumberOfIds() / (nextCell.dim+1);
                            for(int i=0;i<(nextCell.ntri*(nextCell.dim+1));i++)
                            {
                                vtkIdType j = ptIds->GetId(i); // cell ids have been set with local ids
                                DBG_ASSERT( j>=0 && j<nextCell.np );
                                nextCell.triangulation[i] = j;
                            }
                        }
                        else
                        {
                            debug3<<"Triangulation failed\n";
                        }
                        nextCell.needTriangulation = false;
                        vtkcell = 0;
                    }

                    // switch to next cell
                    if (!twoMaterialOptimization)
                    {
                        cell = nextCell;
                    }

                } // end of 'interface was found'
                else
                {
                    vtkcell = 0;
                }

            } // end of 'cell is ok'

            else // cell is ignored
            {
                //vtkWarningMacro(<<"ignoring cell #"<<ci<<", m="<<m<<", mi="<<mi<<", frac="<<fraction<<"\n");
            }

            // update reference volume
            if( !twoMaterialOptimization )
            {
                referenceVolume -= fraction;
            }

        } // for materials

    } // for cells
    delete [] pointArrayOffset;
    delete [] inPointArrays;
    delete [] inCellArrays;

    ptIds->Delete();
    pts->Delete();
    cpsCell->Delete();
    delete [] interpolatedValues;
    delete [] matOrdering;

    // finish output creation
    for(int m=0;m<nmat;m++)
    {
        if( Mats[m].cellCount > Mats[m].numberOfCells )
        {
            debug3<<"Mat #"<<m<<" : cellCount="<<Mats[m].cellCount<<", numberOfCells="<<Mats[m].numberOfCells<<"\n";
        }
        if( Mats[m].pointCount > Mats[m].numberOfPoints )
        {
            debug3<<"Mat #"<<m<<" : pointCount="<<Mats[m].pointCount<<", numberOfPoints="<<Mats[m].numberOfPoints<<"\n";
        }

        delete [] Mats[m].pointMap;

        vtkUnstructuredGrid* ugOutput = vtkUnstructuredGrid::New();

        // set points
        Mats[m].outPointArrays[nPointData-1]->Squeeze();
        vtkPoints* points = vtkPoints::New(origCoordType);
        points->SetDataTypeToDouble();
        //points->SetDataTypeToFloat();
        points->SetNumberOfPoints( Mats[m].pointCount );
        points->SetData( Mats[m].outPointArrays[nPointData-1] );
        Mats[m].outPointArrays[nPointData-1]->Delete();
        ugOutput->SetPoints( points );
        points->Delete();

        // set cell connectivity
        //
        // Mats[m].cells is arranged in the old VTK (pre v9) way:
        //   c1_nids, c1_id_1, c1_id_2, ... c1_id_n-1, c2_nids, c2_id_1 ...
        //
        // VTK (v9+) connectivity array only contains the point ids for each cell
        // with offsets into that array separately and cellLocations aren't needed

        int conn_size = Mats[m].cellArrayCount - Mats[m].cellCount;

        vtkNew<vtkIdTypeArray> offsets;
        offsets->SetNumberOfValues(Mats[m].cellCount+1);
        vtkIdType oid = 0;
        offsets->SetValue(oid++, 0);

        vtkNew<vtkIdTypeArray> connectivity;
        connectivity->SetNumberOfValues(conn_size);

        const vtkIdType *data = &Mats[m].cells[0];
        const vtkIdType * const dEnd = data + Mats[m].cellArrayCount;
        vtkIdType offset = 0;
        vtkIdType cid = 0;
        while(data < dEnd)
        {
            vtkIdType numPts = *data++;
            offset += numPts;
            offsets->SetValue(oid++, offset);
            while(numPts-- > 0)
            {
                connectivity->SetValue(cid++, *data++);
            }
        }

        vtkCellArray* cellArray = vtkCellArray::New();
        cellArray->SetData(offsets, connectivity);

        // set cell types
        vtkUnsignedCharArray *cellTypes = vtkUnsignedCharArray::New();
        cellTypes->SetNumberOfValues( Mats[m].cellCount );
        unsigned char* cellTypesPtr = cellTypes->WritePointer(0,Mats[m].cellCount);
        for(vtkIdType i=0;i<Mats[m].cellCount;i++)
            cellTypesPtr[i] = Mats[m].cellTypes[i];

        // attach conectivity arrays to data set
        ugOutput->SetCells( cellTypes, cellArray );
        cellArray->Delete();
        cellTypes->Delete();

        // attach point arrays
        for(int i=0;i<nPointData-1;i++)
        {
            Mats[m].outPointArrays[i]->Squeeze();
            ugOutput->GetPointData()->AddArray( Mats[m].outPointArrays[i] );
            Mats[m].outPointArrays[i]->Delete();
        }

        // attach cell arrays
        for(int i=0;i<nCellData;i++)
        {
            Mats[m].outCellArrays[i]->Squeeze();
            ugOutput->GetCellData()->AddArray( Mats[m].outCellArrays[i] );
            Mats[m].outCellArrays[i]->Delete();
        }

        delete [] Mats[m].outCellArrays;
        delete [] Mats[m].outPointArrays;

        // activate attributes similarily to input
        for(int i=0;i<vtkDataSetAttributes::NUM_ATTRIBUTES;i++)
        {
            vtkDataArray* attr = input->GetCellData()->GetAttribute(i);
            if( attr!=0 )
            {
                ugOutput->GetCellData()->SetActiveAttribute(attr->GetName(),i);
            }
        }
        for(int i=0;i<vtkDataSetAttributes::NUM_ATTRIBUTES;i++)
        {
            vtkDataArray* attr = input->GetPointData()->GetAttribute(i);
            if( attr!=0 )
            {
                ugOutput->GetPointData()->SetActiveAttribute(attr->GetName(),i);
            }
        }

        // add material data set to multiblock output
        output[m] = ugOutput;
        //ugOutput->Delete();
    }
    delete [] Mats;

    return 1;
}


#undef GET_POINT_DATA

